#ifndef POISSONOP_H
#define POISSONOP_H

#include "Core/RowSparse.h"
#include "Embedded/EmbeddedGrid.h"

template <int Dim>
class PoissonOp
{
public:
  using iVec = Vec<int, Dim>;
  using rVec = Vec<Real, Dim>;
  using EmGrid = EmbeddedGrid<Dim>;
  enum { fitOrder = 5 };

  struct LaplacianDescriptor;

  PoissonOp(const EmGrid &aGrid, const char *aBcTypes);

  void apply(const Tensor<Real, Dim> &aData, Tensor<Real, Dim> &outData) const;

  void apply(const Tensor<Real, Dim> &aData, const Tensor<Real, Dim> &rhs, Tensor<Real, Dim> &rsd) const;

public:
  const EmGrid &emGrid;
  char  bcTypes[2*Dim+2];
  RowSparse<iVec> matOp;

  const Real dangerousCN = 2e4;
};

//======================================================================

template <int Dim>
class PoissonBottomSolver
{
public:
  using iVec = Vec<int, Dim>;

  PoissonBottomSolver(const PoissonOp<Dim> &aPoissonOp);

  void solve(Tensor<Real, Dim> &aPhi, const Tensor<Real, Dim> &rhs, int) const;

//  void solve(Tensor<Real, Dim> &aPhi, const Tensor<Real, Dim> &rhs, const Tensor<Real, Dim-1> *bcData) const;

protected:
  template <class T>
  void doSolve(Tensor<Real, Dim> &aPhi, const Tensor<Real, Dim> &rhs, const T &bcData) const;

protected:
  const EmbeddedGrid<Dim> &emGrid;
  const PoissonOp<Dim> &poissonOp;
  bool            pureNeumann;
  int             numUnknowns;
  int             numActualUnknowns;
  const char *bcTypes;
  Tensor<int, Dim> unknownMap;
  Tensor<Real, 2> poissonMat;
  Tensor<Real, 2> poissonFact;
  Tensor<int, 1>  ipiv;
};

#endif //POISSONOP_H
